clear
discard
set more off

adopath ++ "./ado/"
adopath ++ "../tools/"
adopath ++ "../ssc/"

******************
*** Parameters ***
******************

local file_name "./input/individual_dataset.csv" // Data file
local output_folder "./output/idiosyncratic" // Output directory
local K = 4 // Number of lags to include when calculating local projections
local L = 16 // Number of IRF coefficients to estimate
local alpha = 0.05 // Significance level

*************************************
*** Import data and set data type ***
*************************************

import delimited using `file_name', case(preserve)

gen qtr = quarterly(DATE,"YQ")

**************************
*** Generate variables ***
**************************

* Forecast errors
gen Efor_Step2 = Realiz1 - SPFfor_Step2

* Deviations from consensus
bysort DATE: egen Consensusfor_Step2 = mean(SPFfor_Step2)
gen DevFromConsensusfor_Step2 = SPFfor_Step2 - Consensusfor_Step2

* Set to panel data
xtset ID qtr

* Generate winsorized variables
winsor2 Efor_Step2 DevFromConsensusfor_Step2 SPFfor_Step2, suffix(_w) cuts(5 95)

***********************************
*** Calculate local projections ***
***********************************

* Estimate a simple AR(4) for deviations from consensus
reghdfe DevFromConsensusfor_Step2 l1.DevFromConsensusfor_Step2 l2.DevFromConsensusfor_Step2 l3.DevFromConsensusfor_Step2 l4.DevFromConsensusfor_Step2, /// 
		absorb(ID) cluster(ID qtr)
matrix AR_coefs = e(b)
matrix AR_coefs_cov = e(V)

* Estimate with deviations from consensus as IV
matrix res_baseline = J(`L', 8, .)
matrix colnames res_baseline = lag IRF se lb ub df_resid K F_stat
	
forvalues ii = 1/`L' {
	local ii_1 = `ii' - 1
	quietly ivreghdfe f`ii_1'.Efor_Step2 l1.DevFromConsensusfor_Step2 l2.DevFromConsensusfor_Step2 l3.DevFromConsensusfor_Step2 l4.DevFromConsensusfor_Step2 /// 
		(SPFfor_Step2 = DevFromConsensusfor_Step2), /// 
		absorb(ID) cluster(ID qtr)
		
	* Collect the results and calculate 
	* confidence intervals at the given significance level
	local se = _se[SPFfor_Step2]
	local df_r = e(df_r) 
	local t_crit = invt(`df_r', `alpha' / 2)
	
	matrix res_baseline [`ii', 1] = `ii'
	matrix res_baseline [`ii', 2] = _b[SPFfor_Step2]
	matrix res_baseline [`ii', 3] = `se'
	matrix res_baseline [`ii', 4] = _b[SPFfor_Step2] + ///
		`t_crit' * `se'
	matrix res_baseline [`ii', 5] = _b[SPFfor_Step2] - ///
		`t_crit' * `se'
	matrix res_baseline [`ii', 6] = `df_r'
	matrix res_baseline [`ii', 7] = `K'
	matrix res_baseline [`ii', 8] = e(widstat)
}	

* Estimate with deviations from consensus as IV -- winsorized
matrix res_winsorized = J(`L', 8, .)
matrix colnames res_winsorized = lag IRF se lb ub df_resid K F_stat
	
forvalues ii = 1/`L' {
	local ii_1 = `ii' - 1
	quietly ivreghdfe f`ii_1'.Efor_Step2_w l1.DevFromConsensusfor_Step2_w l2.DevFromConsensusfor_Step2_w l3.DevFromConsensusfor_Step2_w l4.DevFromConsensusfor_Step2_w /// 
		(SPFfor_Step2_w = DevFromConsensusfor_Step2_w), /// 
		absorb(ID) cluster(ID qtr)
		
	* Collect the results and calculate 
	* confidence intervals at the given significance level
	local se = _se[SPFfor_Step2_w]
	local df_r = e(df_r) 
	local t_crit = invt(`df_r', `alpha' / 2)
	
	matrix res_winsorized [`ii', 1] = `ii'
	matrix res_winsorized [`ii', 2] = _b[SPFfor_Step2_w]
	matrix res_winsorized [`ii', 3] = `se'
	matrix res_winsorized [`ii', 4] = _b[SPFfor_Step2_w] + ///
		`t_crit' * `se'
	matrix res_winsorized [`ii', 5] = _b[SPFfor_Step2_w] - ///
		`t_crit' * `se'
	matrix res_winsorized [`ii', 6] = `df_r'
	matrix res_winsorized [`ii', 7] = `K'
	matrix res_winsorized [`ii', 8] = e(widstat)
}	

* Estimate with JK shocks as IV
matrix res_JK_shocks = J(`L', 8, .)
matrix colnames res_JK_shocks = lag IRF se lb ub df_resid K F_stat
	
forvalues ii = 1/`L' {
	local ii_1 = `ii' - 1
	quietly ivreghdfe f`ii_1'.Efor_Step2 l1.JK_shock l2.JK_shock l3.JK_shock l4.JK_shock /// 
		(SPFfor_Step2 = JK_shock), /// 
		absorb(ID) cluster(ID qtr)
		
	* Collect the results and calculate 
	* confidence intervals at the given significance level
	local se = _se[SPFfor_Step2]
	local df_r = e(df_r) 
	local t_crit = invt(`df_r', `alpha' / 2)
	
	matrix res_JK_shocks [`ii', 1] = `ii'
	matrix res_JK_shocks [`ii', 2] = _b[SPFfor_Step2]
	matrix res_JK_shocks [`ii', 3] = `se'
	matrix res_JK_shocks [`ii', 4] = _b[SPFfor_Step2] + ///
		`t_crit' * `se'
	matrix res_JK_shocks [`ii', 5] = _b[SPFfor_Step2] - ///
		`t_crit' * `se'
	matrix res_JK_shocks [`ii', 6] = `df_r'
	matrix res_JK_shocks [`ii', 7] = `K'
	matrix res_JK_shocks [`ii', 8] = e(widstat)
}

******************
*** Get graphs ***
******************

* Clean up
shell rm -r `output_folder'
shell mkdir -p ./output
shell mkdir `output_folder'

* Get plots
plot_IRF "res_baseline" `output_folder' "res_baseline"
plot_IRF "res_JK_shocks" `output_folder' "res_JK_shocks"

* Save estimated IRFs
mat2txt, matrix(res_baseline) saving(`output_folder'/IRF_idiosyncratic.txt)
mat2txt, matrix(res_winsorized) saving(`output_folder'/IRF_idiosyncratic_winsorized.txt)
mat2txt, matrix(res_JK_shocks) saving(`output_folder'/IRF_idiosyncratic_JK_shocks.txt)
mat2txt, matrix(AR_coefs) saving(`output_folder'/AR_coefs_dev_consensus.txt)
mat2txt, matrix(AR_coefs_cov) saving(`output_folder'/AR_coefs_cov_dev_consensus.txt)
